options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -3212.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -44.2348   0.000571428   0.000101815           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -44.23
##    b0         6.31
##    b1         1.87
##    s          5.54
##    m0[1]     18.23
##    m0[2]     21.47
##    m0[3]     13.17
##    m0[4]     22.25
##    m0[5]      6.33
##    m0[6]     10.94
##    m0[7]     43.59
##    m0[8]     16.32
##    m0[9]      5.80
##    m0[10]    33.11
##    m0[11]    26.81
##    m0[12]    41.24
##    m0[13]    17.51
##    m0[14]    28.43
##    m0[15]    28.94
##    m0[16]    40.29
##    m0[17]    24.08
##    m0[18]    21.39
##    m0[19]    15.65
##    m0[20]    40.46
##    y0[1]     22.51
##    y0[2]     29.22
##    y0[3]     14.21
##    y0[4]     22.20
##    y0[5]     11.34
##    y0[6]      6.84
##    y0[7]     46.17
##    y0[8]     17.39
##    y0[9]     10.52
##    y0[10]    36.19
##    y0[11]    28.28
##    y0[12]    50.55
##    y0[13]    25.44
##    y0[14]    20.28
##    y0[15]    30.99
##    y0[16]    29.41
##    y0[17]    21.73
##    y0[18]    13.53
##    y0[19]     6.07
##    y0[20]    45.62
##    m1[1]      5.80
##    m1[2]     10.00
##    m1[3]     14.20
##    m1[4]     18.39
##    m1[5]     22.59
##    m1[6]     26.79
##    m1[7]     30.99
##    m1[8]     35.19
##    m1[9]     39.39
##    m1[10]    43.59
##    y1[1]      8.46
##    y1[2]     13.16
##    y1[3]     15.73
##    y1[4]     15.14
##    y1[5]     19.36
##    y1[6]     30.33
##    y1[7]     26.02
##    y1[8]     40.43
##    y1[9]     43.25
##    y1[10]    47.33
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -44.18 -43.76 1.44 1.08 -47.10 -42.70 1.00      405      425
##    b0       6.70   6.56 2.71 2.19   2.40  11.36 1.01      265      204
##    b1       1.83   1.83 0.24 0.22   1.43   2.24 1.01      357      409
##    s        6.32   6.16 1.20 1.05   4.69   8.53 1.00      951      802
##    m0[1]   18.39  18.40 1.65 1.45  15.68  21.24 1.01      336      335
##    m0[2]   21.58  21.62 1.50 1.41  19.15  24.00 1.01      458      518
##    m0[3]   13.43  13.38 2.03 1.73  10.16  16.97 1.01      276      225
##    m0[4]   22.34  22.39 1.49 1.40  20.00  24.73 1.01      511      590
##    m0[5]    6.72   6.58 2.71 2.19   2.42  11.39 1.01      265      204
##    m0[6]   11.24  11.15 2.24 1.87   7.62  15.08 1.01      269      217
##    m0[7]   43.27  43.27 2.88 2.80  38.51  47.97 1.00      999      933
##    m0[8]   16.52  16.50 1.78 1.52  13.72  19.58 1.01      303      256
##    m0[9]    6.20   6.06 2.77 2.24   1.82  10.95 1.01      264      204
##    m0[10]  33.00  32.98 1.86 1.80  29.88  36.00 1.00     2059     1320
##    m0[11]  26.81  26.83 1.50 1.45  24.30  29.24 1.00     1240     1072
##    m0[12]  40.97  40.95 2.63 2.55  36.59  45.27 1.00     1190     1114
##    m0[13]  17.69  17.68 1.69 1.47  14.93  20.66 1.01      321      298
##    m0[14]  28.40  28.41 1.56 1.51  25.77  30.90 1.00     1780     1190
##    m0[15]  28.91  28.92 1.59 1.52  26.24  31.44 1.00     1885     1207
##    m0[16]  40.03  40.02 2.53 2.47  35.81  44.16 1.00     1289     1058
##    m0[17]  24.13  24.15 1.46 1.41  21.75  26.54 1.00      698      754
##    m0[18]  21.50  21.54 1.51 1.40  19.06  23.93 1.01      452      518
##    m0[19]  15.87  15.84 1.83 1.56  12.99  19.01 1.01      295      244
##    m0[20]  40.20  40.18 2.55 2.48  35.96  44.37 1.00     1269     1058
##    y0[1]   18.28  18.31 6.72 6.38   7.25  29.63 1.00     1367     1751
##    y0[2]   21.51  21.33 6.60 6.36  10.71  32.22 1.00     2103     1946
##    y0[3]   13.47  13.66 6.72 6.38   2.69  24.19 1.00     1377     1791
##    y0[4]   22.24  22.16 6.65 6.52  11.69  33.54 1.00     1465     1581
##    y0[5]    6.69   6.83 6.96 6.82  -4.66  18.16 1.01      976     1169
##    y0[6]   11.17  11.31 6.82 6.51   0.23  22.15 1.00     1277      934
##    y0[7]   43.55  43.44 7.18 6.42  32.38  55.60 1.00     1714     1791
##    y0[8]   16.29  16.36 6.66 6.23   5.28  27.30 1.00     1279     1661
##    y0[9]    6.21   6.17 6.94 6.47  -5.33  17.54 1.00      875      994
##    y0[10]  32.90  32.88 6.76 6.53  21.66  44.03 1.00     1898     1930
##    y0[11]  26.87  26.91 6.71 6.57  16.33  37.58 1.00     2004     1974
##    y0[12]  41.23  41.02 6.99 6.63  30.25  52.79 1.00     1955     1856
##    y0[13]  17.60  17.79 6.86 6.67   6.45  28.43 1.01     1222     1330
##    y0[14]  28.49  28.47 6.68 6.32  17.66  39.27 1.00     2124     1970
##    y0[15]  28.80  28.59 6.75 6.23  17.62  40.05 1.00     2166     1907
##    y0[16]  40.05  39.94 7.04 6.85  28.60  51.76 1.00     2077     1955
##    y0[17]  24.04  24.01 6.58 6.21  13.28  34.51 1.00     1736     1863
##    y0[18]  21.63  21.61 6.68 6.49  10.88  32.62 1.00     1783     1876
##    y0[19]  15.79  16.00 6.59 6.08   4.72  26.17 1.00     1529     1731
##    y0[20]  40.22  40.23 6.88 6.46  28.71  51.18 1.00     1985     1838
##    m1[1]    6.20   6.06 2.77 2.24   1.82  10.95 1.01      264      204
##    m1[2]   10.32  10.23 2.33 1.93   6.60  14.33 1.01      267      215
##    m1[3]   14.44  14.41 1.94 1.65  11.31  17.85 1.01      282      236
##    m1[4]   18.56  18.56 1.64 1.45  15.86  21.37 1.01      340      335
##    m1[5]   22.68  22.72 1.48 1.40  20.35  25.08 1.01      539      612
##    m1[6]   26.79  26.82 1.50 1.45  24.29  29.23 1.00     1235     1072
##    m1[7]   30.91  30.91 1.71 1.63  28.06  33.71 1.00     2088     1297
##    m1[8]   35.03  35.02 2.03 1.96  31.62  38.36 1.00     1845     1286
##    m1[9]   39.15  39.14 2.44 2.36  35.07  43.10 1.00     1402     1037
##    m1[10]  43.27  43.27 2.88 2.80  38.51  47.97 1.00      999      933
##    y1[1]    6.25   6.31 6.84 6.60  -5.03  17.77 1.00     1070     1388
##    y1[2]   10.18  10.37 6.81 6.66  -1.02  21.67 1.00     1019     1279
##    y1[3]   14.50  14.47 6.77 6.52   3.63  25.58 1.00     1171     1522
##    y1[4]   18.65  18.57 6.82 6.49   7.73  29.42 1.00     1836     1923
##    y1[5]   22.57  22.71 6.85 6.44  11.53  33.82 1.00     1863     1933
##    y1[6]   26.86  26.92 6.49 6.39  16.25  37.07 1.00     2135     2009
##    y1[7]   30.89  30.80 6.72 6.35  20.04  41.83 1.00     1600     1907
##    y1[8]   34.92  34.99 6.91 6.41  23.67  46.13 1.00     1970     1894
##    y1[9]   39.28  39.23 6.77 6.34  27.99  50.30 1.00     1933     1953
##    y1[10]  43.37  43.53 7.14 6.89  31.65  54.73 1.00     1618     1448

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.6 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 finished in 0.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 0.7 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -57.13 -58.46 10.96 9.52 -73.12 -36.69 1.02      147      159
##    b0       6.73   6.67  2.86 2.77   2.26  11.46 1.01      873     1170
##    b1       1.82   1.82  0.26 0.24   1.40   2.24 1.01      728      979
##    s        4.36   4.42  1.93 2.11   1.24   7.49 1.06       65      171
##    sx       2.36   2.37  1.27 1.29   0.42   4.47 1.04      109      148
##    x0[1]    4.98   5.22  1.87 1.83   1.87   7.59 1.02      326     1386
##    x0[2]    9.42   9.32  1.72 1.76   6.82  12.30 1.02      384     1902
##    x0[3]    1.94   2.28  2.11 2.06  -1.66   4.75 1.02      286      883
##    x0[4]    8.74   8.69  1.54 1.19   6.21  11.34 1.01     3254     1758
##    x0[5]    1.00   0.93  1.80 1.67  -1.76   3.89 1.01      846     1874
##    x0[6]    2.63   2.63  1.67 1.38  -0.14   5.31 1.01     1881     1429
##    x0[7]   18.38  18.54  1.96 2.08  15.11  21.20 1.01      511     1670
##    x0[8]    5.62   5.59  1.54 1.26   3.22   8.04 1.00     2340     1800
##    x0[9]   -0.59  -0.42  1.75 1.41  -3.71   2.05 1.00     1395      907
##    x0[10]  14.19  14.20  1.56 1.33  11.68  16.82 1.01     2232     1896
##    x0[11]  14.01  13.83  2.63 3.18  10.51  18.30 1.04      146     1185
##    x0[12]  18.40  18.43  1.66 1.42  15.77  21.07 1.00     1382     1201
##    x0[13]   1.81   2.13  3.39 4.15  -3.84   6.26 1.04      121      692
##    x0[14]  13.38  13.17  1.89 1.91  10.70  16.65 1.02      247     1281
##    x0[15]  12.49  12.40  1.54 1.25  10.10  15.12 1.00     1669     1386
##    x0[16]  18.49  18.30  1.71 1.40  16.00  21.46 1.01     1338     1488
##    x0[17]  10.32  10.20  1.64 1.42   7.90  13.01 1.01      900     2024
##    x0[18]   9.92   9.83  1.97 2.20   7.05  13.16 1.03      146     1397
##    x0[19]   3.62   3.83  1.88 1.79   0.38   6.25 1.02      305     1202
##    x0[20]  18.74  18.54  1.85 1.41  16.07  21.85 1.01     1256     1069
##    m0[1]   15.93  15.80  3.16 3.29  11.08  21.04 1.01      507     2554
##    m0[2]   23.85  23.88  3.09 3.27  18.84  28.70 1.02      338     2115
##    m0[3]   10.43  10.27  3.56 3.70   5.02  16.32 1.02      377     2073
##    m0[4]   22.64  22.66  2.72 2.33  18.14  27.02 1.01     4145     2190
##    m0[5]    8.68   8.88  3.38 3.41   3.03  13.77 1.01      585     2142
##    m0[6]   11.64  11.57  2.96 2.62   6.84  16.56 1.01     2756     1891
##    m0[7]   40.06  39.83  3.74 3.92  34.63  46.38 1.02      291     1585
##    m0[8]   17.02  17.05  2.75 2.42  12.57  21.48 1.00     2800     2598
##    m0[9]    5.84   5.71  3.11 2.74   0.86  11.10 1.00     2283     1892
##    m0[10]  32.49  32.45  2.81 2.46  27.93  37.22 1.00     3554     2242
##    m0[11]  32.09  31.99  4.50 5.66  25.17  38.85 1.03      113     1931
##    m0[12]  40.07  40.01  2.99 2.65  35.18  45.05 1.00     1719     1999
##    m0[13]  10.23  10.22  5.74 7.51   1.82  18.94 1.04      128     1210
##    m0[14]  30.98  31.05  3.22 3.50  25.76  36.05 1.02      308     2435
##    m0[15]  29.41  29.46  2.72 2.43  24.92  33.67 1.00     3499     2413
##    m0[16]  40.22  40.25  2.94 2.73  35.33  45.08 1.00     3383     2444
##    m0[17]  25.48  25.56  2.88 2.70  20.81  30.07 1.01     1079     2531
##    m0[18]  24.75  24.75  3.52 3.88  19.20  30.23 1.03      159     2588
##    m0[19]  13.45  13.28  3.23 3.22   8.46  18.74 1.01      546     2091
##    m0[20]  40.66  40.75  3.13 2.79  35.26  45.71 1.00     3493     2437
##    y0[1]   15.85  15.14  5.75 4.86   7.33  26.20 1.01     1876     2915
##    y0[2]   24.02  24.53  5.71 4.99  13.90  32.78 1.01     1560     2453
##    y0[3]   10.44   9.68  5.81 5.21   2.02  20.76 1.01     1236     3227
##    y0[4]   22.49  22.63  5.44 4.61  13.40  31.26 1.01     3869     2910
##    y0[5]    8.64   9.08  5.99 5.07  -2.08  17.55 1.01     1567     2687
##    y0[6]   11.62  11.63  5.60 4.74   2.19  20.64 1.01     4354     3407
##    y0[7]   39.93  39.16  6.10 5.39  31.26  50.63 1.01     1177     1901
##    y0[8]   17.09  17.15  5.39 4.49   8.10  26.07 1.01     3307     3332
##    y0[9]    5.84   5.74  5.66 4.60  -3.47  15.43 1.01     3313     2505
##    y0[10]  32.59  32.44  5.53 4.48  23.79  42.04 1.01     3861     3171
##    y0[11]  31.92  32.80  6.55 6.44  20.11  40.84 1.02      289     2071
##    y0[12]  40.09  39.79  5.60 4.51  30.91  49.79 1.01     2897     3274
##    y0[13]  10.25   9.20  7.55 7.87   0.36  23.75 1.02      264     1353
##    y0[14]  30.84  31.65  5.78 4.94  20.55  39.18 1.01     1273     2721
##    y0[15]  29.39  29.58  5.48 4.44  20.04  38.04 1.01     4142     2891
##    y0[16]  40.18  40.20  5.51 4.63  31.23  49.36 1.01     3823     3145
##    y0[17]  25.60  26.03  5.52 4.66  15.94  33.97 1.01     2616     2582
##    y0[18]  24.65  25.49  5.96 5.31  13.92  33.43 1.01      986     2472
##    y0[19]  13.35  12.70  5.74 4.82   4.92  23.58 1.01     1807     2595
##    y0[20]  40.62  40.91  5.80 4.79  30.85  49.71 1.01     3553     2938
##    m1[1]    6.23   6.17  2.92 2.82   1.67  11.07 1.01      866     1135
##    m1[2]   10.33  10.29  2.44 2.33   6.44  14.34 1.01      936     1265
##    m1[3]   14.42  14.38  2.01 1.92  11.15  17.79 1.00     1080     1762
##    m1[4]   18.52  18.51  1.66 1.58  15.79  21.23 1.00     1356     2225
##    m1[5]   22.61  22.61  1.47 1.39  20.18  24.96 1.00     1715     2550
##    m1[6]   26.70  26.70  1.50 1.44  24.26  29.08 1.00     1658     2398
##    m1[7]   30.80  30.78  1.72 1.66  27.93  33.56 1.00     1229     1923
##    m1[8]   34.89  34.87  2.09 2.04  31.47  38.24 1.01      975     1417
##    m1[9]   38.99  38.96  2.54 2.44  34.84  43.11 1.01      866     1304
##    m1[10]  43.08  43.03  3.02 2.87  38.15  48.07 1.01      814     1295
##    x1[1]   -0.27  -0.26  2.71 1.89  -4.68   4.12 1.01     3822     3082
##    x1[2]    2.01   1.98  2.65 1.82  -2.10   6.32 1.02     3924     2630
##    x1[3]    4.23   4.24  2.66 1.91  -0.24   8.72 1.02     4081     1948
##    x1[4]    6.42   6.46  2.65 1.88   2.01  10.71 1.01     3486     2368
##    x1[5]    8.71   8.70  2.59 1.83   4.46  13.06 1.01     4154     2681
##    x1[6]   10.94  10.97  2.72 1.95   6.42  15.17 1.02     4062     2458
##    x1[7]   13.23  13.22  2.63 1.96   8.91  17.63 1.01     4086     2501
##    x1[8]   15.45  15.46  2.66 1.94  11.05  19.77 1.01     3940     2459
##    x1[9]   17.71  17.72  2.71 1.95  13.39  22.06 1.01     4017     2355
##    x1[10]  19.94  19.94  2.65 1.86  15.59  24.18 1.01     4498     2475
##    y1[1]    6.24   6.24  5.67 4.79  -2.96  15.20 1.01     1932     2447
##    y1[2]   10.46  10.49  5.25 4.46   1.83  19.13 1.01     2863     2619
##    y1[3]   14.48  14.51  5.14 4.31   6.07  23.02 1.01     3040     2695
##    y1[4]   18.38  18.52  5.13 4.20   9.80  26.61 1.01     2842     2651
##    y1[5]   22.73  22.74  5.00 4.04  14.65  31.07 1.01     3721     3043
##    y1[6]   26.60  26.59  4.99 3.98  18.45  34.79 1.01     3139     2788
##    y1[7]   30.70  30.70  4.96 4.22  22.34  38.87 1.01     3151     2833
##    y1[8]   34.96  34.86  5.32 4.48  26.42  44.01 1.01     2489     2777
##    y1[9]   38.99  38.97  5.26 4.47  30.32  47.55 1.01     2444     3103
##    y1[10]  43.16  43.23  5.72 5.00  33.58  52.59 1.01     1520     2616

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -2635.41 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -33.4085   0.000232464   0.000110767      0.3378      0.9899       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -33.41
##    b0         6.76
##    b1         1.80
##    s          3.22
##    m0[1]     12.15
##    m0[2]     22.66
##    m0[3]     16.16
##    m0[4]     17.89
##    m0[5]     16.68
##    m0[6]     21.87
##    m0[7]     22.46
##    m0[8]     18.49
##    m0[9]     15.72
##    m0[10]     8.08
##    m0[11]    15.38
##    m0[12]    17.28
##    m0[13]    20.49
##    m0[14]     9.05
##    m0[15]    12.83
##    m0[16]    15.95
##    m0[17]    19.81
##    m0[18]    21.84
##    m0[19]     9.24
##    m0[20]    21.01
##    y0[1]     12.49
##    y0[2]     21.55
##    y0[3]     16.14
##    y0[4]     13.73
##    y0[5]     18.04
##    y0[6]     23.37
##    y0[7]     22.50
##    y0[8]     19.11
##    y0[9]     17.11
##    y0[10]     8.63
##    y0[11]    11.79
##    y0[12]    18.73
##    y0[13]    17.70
##    y0[14]     5.52
##    y0[15]    13.27
##    y0[16]    15.39
##    y0[17]    12.43
##    y0[18]    18.51
##    y0[19]    11.25
##    y0[20]    17.72
##    m1[1]      8.08
##    m1[2]      9.70
##    m1[3]     11.32
##    m1[4]     12.94
##    m1[5]     14.56
##    m1[6]     16.18
##    m1[7]     17.80
##    m1[8]     19.42
##    m1[9]     21.04
##    m1[10]    22.66
##    y1[1]      7.39
##    y1[2]      9.52
##    y1[3]     14.49
##    y1[4]     15.14
##    y1[5]     12.61
##    y1[6]     18.93
##    y1[7]      8.08
##    y1[8]     22.46
##    y1[9]     22.23
##    y1[10]    20.24
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.91 -33.59 1.33 1.14 -36.51 -32.41 1.00      550      777
##    b0       6.70   6.70 2.09 2.06   3.27   9.97 1.01      253      326
##    b1       1.81   1.80 0.34 0.34   1.25   2.36 1.01      300      418
##    s        3.69   3.58 0.67 0.60   2.77   4.94 1.00     1290     1290
##    m0[1]   12.12  12.11 1.23 1.20  10.13  14.11 1.01      286      379
##    m0[2]   22.69  22.69 1.40 1.35  20.43  25.04 1.00     1044     1332
##    m0[3]   16.16  16.17 0.87 0.85  14.74  17.62 1.00      801     1018
##    m0[4]   17.89  17.92 0.89 0.86  16.44  19.34 1.00     1869     1569
##    m0[5]   16.68  16.70 0.86 0.84  15.24  18.12 1.00     1033     1350
##    m0[6]   21.89  21.90 1.29 1.24  19.82  24.04 1.00     1256     1423
##    m0[7]   22.49  22.50 1.37 1.32  20.28  24.81 1.00     1088     1332
##    m0[8]   18.50  18.52 0.92 0.90  17.02  20.01 1.00     2244     1623
##    m0[9]   15.71  15.72 0.89 0.86  14.27  17.18 1.01      652      987
##    m0[10]   8.03   8.06 1.86 1.83   4.99  10.91 1.01      253      329
##    m0[11]  15.37  15.38 0.90 0.87  13.92  16.91 1.01      563      817
##    m0[12]  17.28  17.30 0.87 0.83  15.82  18.70 1.00     1396     1415
##    m0[13]  20.51  20.52 1.11 1.07  18.69  22.32 1.00     1838     1545
##    m0[14]   9.00   9.02 1.70 1.67   6.21  11.65 1.01      255      314
##    m0[15]  12.81  12.78 1.14 1.13  10.96  14.68 1.01      307      511
##    m0[16]  15.94  15.95 0.88 0.85  14.52  17.41 1.00      729     1029
##    m0[17]  19.83  19.85 1.03 0.99  18.15  21.50 1.00     2124     1556
##    m0[18]  21.87  21.87 1.29 1.24  19.80  24.01 1.00     1263     1444
##    m0[19]   9.20   9.22 1.67 1.63   6.45  11.79 1.01      256      314
##    m0[20]  21.03  21.03 1.17 1.13  19.11  22.97 1.00     1609     1514
##    y0[1]   12.10  11.99 4.01 3.85   5.63  18.87 1.00     1336     1448
##    y0[2]   22.59  22.71 4.03 3.91  15.84  28.89 1.00     1987     1845
##    y0[3]   16.10  16.08 3.82 3.70   9.92  22.39 1.00     1727     1852
##    y0[4]   17.86  17.90 3.97 3.87  11.46  24.29 1.00     2002     1995
##    y0[5]   16.75  16.69 3.92 3.69  10.61  23.27 1.00     1773     1895
##    y0[6]   21.98  21.95 4.17 3.95  15.24  28.79 1.00     1692     1592
##    y0[7]   22.53  22.48 3.95 3.90  16.17  29.34 1.00     1917     1778
##    y0[8]   18.53  18.53 3.86 3.79  12.17  24.77 1.00     1918     2053
##    y0[9]   15.65  15.52 3.82 3.73   9.44  21.88 1.00     2156     1713
##    y0[10]   8.14   8.13 4.17 4.04   1.39  15.20 1.00     1115     1330
##    y0[11]  15.17  15.24 3.81 3.73   8.86  21.35 1.00     1730     1972
##    y0[12]  17.46  17.44 3.85 3.72  11.30  23.98 1.00     1990     1912
##    y0[13]  20.54  20.50 3.94 3.73  14.35  27.27 1.00     2021     1725
##    y0[14]   8.94   8.99 3.99 4.02   2.36  15.25 1.00      835     1614
##    y0[15]  12.78  12.81 3.95 3.79   6.30  18.98 1.00     1243     1744
##    y0[16]  15.84  15.88 3.71 3.59   9.89  22.04 1.00     1991     1681
##    y0[17]  20.01  20.00 3.88 3.70  13.65  26.25 1.00     2029     1916
##    y0[18]  21.94  21.92 3.94 3.72  15.82  28.36 1.00     2009     1934
##    y0[19]   9.19   9.18 4.07 4.02   2.36  15.87 1.00      894     1333
##    y0[20]  20.99  20.85 4.00 3.85  14.76  27.58 1.00     1831     1940
##    m1[1]    8.03   8.06 1.86 1.83   4.99  10.91 1.01      253      329
##    m1[2]    9.66   9.68 1.60 1.56   7.04  12.15 1.01      258      341
##    m1[3]   11.29  11.28 1.35 1.32   9.10  13.46 1.01      271      361
##    m1[4]   12.92  12.89 1.13 1.12  11.10  14.77 1.01      312      501
##    m1[5]   14.55  14.54 0.96 0.91  12.97  16.16 1.01      424      640
##    m1[6]   16.17  16.19 0.87 0.85  14.76  17.64 1.00      807     1075
##    m1[7]   17.80  17.82 0.88 0.85  16.35  19.25 1.00     1790     1568
##    m1[8]   19.43  19.45 1.00 0.96  17.85  21.06 1.00     2248     1543
##    m1[9]   21.06  21.06 1.18 1.13  19.14  23.00 1.00     1597     1514
##    m1[10]  22.69  22.69 1.40 1.35  20.43  25.04 1.00     1044     1332
##    y1[1]    7.99   8.06 4.15 3.91   1.35  14.76 1.00     1000     1345
##    y1[2]    9.62   9.60 4.07 3.86   3.01  16.46 1.00     1122     1682
##    y1[3]   11.34  11.31 4.02 3.99   4.81  17.89 1.00     1308     1694
##    y1[4]   12.93  12.83 4.03 3.75   6.52  19.50 1.00     1607     1841
##    y1[5]   14.47  14.45 3.87 3.76   7.96  20.93 1.00     1727     2012
##    y1[6]   16.21  16.25 3.84 3.56  10.01  22.41 1.00     1884     1800
##    y1[7]   17.68  17.58 3.82 3.79  11.38  24.07 1.00     1964     1754
##    y1[8]   19.45  19.60 3.90 3.77  13.09  25.65 1.00     2213     1971
##    y1[9]   21.07  20.99 3.97 3.62  14.46  27.74 1.00     1983     1775
##    y1[10]  22.67  22.63 4.01 3.91  16.22  29.22 1.00     1972     1971

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -72.2134 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -15.9805   0.000356175    0.00111771           1           1       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -15.98
##    b0         4.64
##    b1         2.10
##    s          0.68
##    m0[1]     10.94
##    m0[2]     23.23
##    m0[3]     15.64
##    m0[4]     17.65
##    m0[5]     16.24
##    m0[6]     22.30
##    m0[7]     22.99
##    m0[8]     18.36
##    m0[9]     15.11
##    m0[10]     6.19
##    m0[11]    14.72
##    m0[12]    16.94
##    m0[13]    20.69
##    m0[14]     7.32
##    m0[15]    11.74
##    m0[16]    15.38
##    m0[17]    19.90
##    m0[18]    22.27
##    m0[19]     7.55
##    m0[20]    21.30
##    y0[1]     11.39
##    y0[2]     22.40
##    y0[3]     16.49
##    y0[4]     17.86
##    y0[5]     -1.41
##    y0[6]     23.63
##    y0[7]     23.74
##    y0[8]     17.73
##    y0[9]     -9.13
##    y0[10]     3.65
##    y0[11]    14.93
##    y0[12]    17.02
##    y0[13]    21.37
##    y0[14]     7.58
##    y0[15]    12.74
##    y0[16]    14.20
##    y0[17]    19.94
##    y0[18]    22.96
##    y0[19]     0.50
##    y0[20]    21.10
##    m1[1]      6.19
##    m1[2]      8.08
##    m1[3]      9.97
##    m1[4]     11.87
##    m1[5]     13.76
##    m1[6]     15.65
##    m1[7]     17.55
##    m1[8]     19.44
##    m1[9]     21.33
##    m1[10]    23.23
##    y1[1]      3.82
##    y1[2]      6.63
##    y1[3]      9.40
##    y1[4]     11.12
##    y1[5]     13.01
##    y1[6]     17.35
##    y1[7]     18.27
##    y1[8]     16.62
##    y1[9]     20.10
##    y1[10]    23.67
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -18.00 -17.60   1.42 1.07 -20.92 -16.51 1.01      589      517
##    b0       4.62   4.61   0.61 0.45   3.63   5.57 1.01      488      366
##    b1       2.10   2.10   0.11 0.10   1.92   2.27 1.01      646      477
##    s        0.85   0.81   0.27 0.27   0.47   1.34 1.00      607      666
##    m0[1]   10.90  10.91   0.35 0.29  10.33  11.44 1.01      583      450
##    m0[2]   23.15  23.16   0.55 0.57  22.24  24.01 1.00     1407     1387
##    m0[3]   15.58  15.60   0.30 0.30  15.05  16.03 1.00     1577      974
##    m0[4]   17.59  17.61   0.34 0.34  17.01  18.12 1.00     2101     1464
##    m0[5]   16.18  16.20   0.31 0.31  15.64  16.65 1.00     1823     1068
##    m0[6]   22.22  22.24   0.51 0.53  21.38  23.02 1.00     1546     1403
##    m0[7]   22.91  22.93   0.54 0.56  22.02  23.76 1.00     1440     1450
##    m0[8]   18.30  18.31   0.36 0.36  17.68  18.85 1.00     2092     1476
##    m0[9]   15.06  15.08   0.30 0.29  14.53  15.50 1.00     1379      857
##    m0[10]   6.16   6.16   0.53 0.40   5.29   7.00 1.01      487      349
##    m0[11]  14.67  14.69   0.30 0.29  14.15  15.12 1.00     1244      747
##    m0[12]  16.89  16.91   0.32 0.32  16.32  17.38 1.00     2026     1527
##    m0[13]  20.62  20.63   0.44 0.46  19.88  21.32 1.00     1832     1564
##    m0[14]   7.29   7.29   0.49 0.37   6.49   8.05 1.01      492      369
##    m0[15]  11.70  11.71   0.33 0.28  11.15  12.19 1.01      639      452
##    m0[16]  15.33  15.35   0.30 0.30  14.81  15.78 1.00     1479      962
##    m0[17]  19.83  19.85   0.41 0.42  19.12  20.47 1.00     1969     1675
##    m0[18]  22.20  22.22   0.51 0.53  21.36  22.99 1.00     1551     1468
##    m0[19]   7.51   7.52   0.48 0.36   6.75   8.26 1.01      494      369
##    m0[20]  21.23  21.24   0.47 0.48  20.44  21.96 1.00     1718     1617
##    y0[1]   11.68  10.89  99.38 1.32   5.91  15.95 1.00     1884     1638
##    y0[2]   22.68  23.17  30.39 1.40  18.04  28.52 1.00     2151     1968
##    y0[3]   15.95  15.59  22.87 1.21  10.46  22.26 1.00     1873     1959
##    y0[4]   17.36  17.62  19.24 1.31  12.63  22.81 1.00     2060     1854
##    y0[5]   16.54  16.16  78.53 1.36  10.38  22.05 1.00     1765     1924
##    y0[6]   30.65  22.25 404.41 1.38  17.01  27.63 1.00     2015     2056
##    y0[7]   22.35  22.95  32.23 1.46  17.76  28.59 1.00     1712     1890
##    y0[8]   18.04  18.27   9.61 1.33  12.62  23.38 1.00     2112     1851
##    y0[9]   17.55  15.07  68.39 1.24  10.32  20.18 1.00     1947     1964
##    y0[10]   5.00   6.17  47.07 1.40   0.78  11.56 1.00     1921     1973
##    y0[11]  14.72  14.64  16.08 1.27   9.16  19.42 1.00     1831     1759
##    y0[12]  16.68  16.90  14.44 1.28  11.67  21.51 1.00     1969     1988
##    y0[13]  20.47  20.61  23.53 1.39  15.16  25.88 1.00     1972     1861
##    y0[14]   2.42   7.34 196.97 1.33   2.66  13.42 1.00     1841     1963
##    y0[15]  13.15  11.66  53.16 1.18   7.51  16.48 1.00     1700     1955
##    y0[16]  15.45  15.37  12.25 1.27  10.70  20.82 1.00     1812     1940
##    y0[17]  19.08  19.90  42.43 1.33  14.72  24.61 1.00     1877     1928
##    y0[18]  23.19  22.19  50.14 1.43  17.12  27.88 1.00     2255     1870
##    y0[19]   7.56   7.54  20.99 1.43   2.53  13.63 1.00     1768     1748
##    y0[20]  18.11  21.30 140.14 1.35  15.61  27.39 1.00     2027     1974
##    m1[1]    6.16   6.16   0.53 0.40   5.29   7.00 1.01      487      349
##    m1[2]    8.05   8.05   0.45 0.35   7.32   8.76 1.01      499      351
##    m1[3]    9.94   9.95   0.38 0.30   9.33  10.52 1.01      537      406
##    m1[4]   11.82  11.84   0.33 0.28  11.28  12.31 1.01      651      456
##    m1[5]   13.71  13.73   0.30 0.27  13.21  14.15 1.01      976      586
##    m1[6]   15.60  15.62   0.30 0.30  15.07  16.05 1.00     1583      959
##    m1[7]   17.49  17.50   0.34 0.34  16.91  18.01 1.00     2140     1465
##    m1[8]   19.37  19.39   0.40 0.41  18.70  19.99 1.00     2041     1645
##    m1[9]   21.26  21.28   0.47 0.48  20.47  22.00 1.00     1712     1621
##    m1[10]  23.15  23.16   0.55 0.57  22.24  24.01 1.00     1407     1387
##    y1[1]    6.82   6.19  55.80 1.37   0.40  11.22 1.00     1635     1644
##    y1[2]    8.06   8.03  30.11 1.26   2.16  13.63 1.00     1891     1767
##    y1[3]    8.79   9.93  29.38 1.31   4.76  15.37 1.00     1936     1890
##    y1[4]   11.18  11.82  22.64 1.25   6.15  16.99 1.00     1856     1940
##    y1[5]   11.07  13.72 181.16 1.25   9.01  18.90 1.00     1791     1713
##    y1[6]   15.33  15.63  34.36 1.27   9.20  21.03 1.00     2103     1872
##    y1[7]   14.07  17.53 187.16 1.28  11.70  22.31 1.00     1840     1882
##    y1[8]   19.89  19.39  25.94 1.30  13.31  25.16 1.00     2141     1846
##    y1[9]   20.55  21.24  32.17 1.36  15.00  26.65 1.00     1929     1763
##    y1[10]  24.10  23.08  42.41 1.44  17.58  28.46 1.00     1973     1972

censored data

y i=1-N, y~N(m,s)  
  acutualy        ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11.stan

objective variable is censored

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -1250.55 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -22.0996   0.000230068   0.000160309           1           1       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -22.10
##    b0        16.93
##    b1         1.47
##    s          3.32
##    m0[1]     23.18
##    m0[2]     22.69
##    m0[3]     18.66
##    m0[4]     23.58
##    m0[5]     23.65
##    m0[6]     19.29
##    m0[7]     26.20
##    m0[8]     21.62
##    m0[9]     28.05
##    m0[10]    28.05
##    m0[11]    26.85
##    m0[12]    26.94
##    m0[13]    27.82
##    y0[1]     25.39
##    y0[2]     23.03
##    y0[3]     18.22
##    y0[4]     21.44
##    y0[5]     27.01
##    y0[6]     16.77
##    y0[7]     22.29
##    y0[8]     20.38
##    y0[9]     27.40
##    y0[10]    28.59
##    y0[11]    32.01
##    y0[12]    32.23
##    y0[13]    26.23
##    m1[1]     17.64
##    m1[2]     18.98
##    m1[3]     20.33
##    m1[4]     21.68
##    m1[5]     23.03
##    m1[6]     24.38
##    m1[7]     25.72
##    m1[8]     27.07
##    m1[9]     28.42
##    m1[10]    29.77
##    y1[1]     17.86
##    y1[2]     21.30
##    y1[3]     24.32
##    y1[4]     21.85
##    y1[5]     19.48
##    y1[6]     16.98
##    y1[7]     24.80
##    y1[8]     24.39
##    y1[9]     31.10
##    y1[10]    36.83
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.70 -22.27 1.65 1.15 -25.85 -21.09 1.02      300      354
##    b0      16.74  16.72 3.33 2.88  11.88  21.37 1.03      263      333
##    b1       1.50   1.51 0.61 0.53   0.60   2.38 1.02      289      390
##    s        4.14   3.94 1.13 0.87   2.80   6.16 1.01      757      491
##    m0[1]   23.11  23.10 1.32 1.19  20.96  25.17 1.01      580      634
##    m0[2]   22.61  22.61 1.40 1.25  20.30  24.76 1.02      464      530
##    m0[3]   18.50  18.47 2.67 2.27  14.50  22.21 1.03      269      331
##    m0[4]   23.51  23.52 1.26 1.14  21.48  25.49 1.01      722      893
##    m0[5]   23.58  23.58 1.26 1.13  21.57  25.55 1.01      753      892
##    m0[6]   19.14  19.13 2.44 2.07  15.44  22.57 1.03      273      334
##    m0[7]   26.18  26.23 1.45 1.29  23.81  28.47 1.00     1399      844
##    m0[8]   21.52  21.53 1.67 1.46  18.88  24.02 1.02      337      326
##    m0[9]   28.06  28.10 1.97 1.74  24.81  31.15 1.01      598      554
##    m0[10]  28.06  28.10 1.97 1.74  24.81  31.15 1.01      599      554
##    m0[11]  26.84  26.90 1.61 1.39  24.23  29.34 1.00      901      704
##    m0[12]  26.94  27.00 1.63 1.43  24.27  29.49 1.00      857      700
##    m0[13]  27.83  27.87 1.90 1.68  24.67  30.77 1.01      633      586
##    y0[1]   23.25  23.24 4.49 4.13  15.98  30.68 1.00     1854     1717
##    y0[2]   22.61  22.69 4.59 4.17  15.13  29.71 1.00     1688     1365
##    y0[3]   18.48  18.60 4.99 4.36  10.57  26.37 1.01      798      885
##    y0[4]   23.47  23.50 4.56 4.22  16.12  30.55 1.00     1836     1799
##    y0[5]   23.62  23.62 4.44 4.16  16.49  30.55 1.00     1767     1702
##    y0[6]   19.26  19.34 4.79 4.48  11.47  26.62 1.01      898     1250
##    y0[7]   26.34  26.32 4.52 4.14  19.10  33.59 1.00     1778     1533
##    y0[8]   21.51  21.62 4.66 4.27  13.41  28.95 1.00     1254     1364
##    y0[9]   27.99  27.96 4.93 4.35  20.15  35.73 1.00     1604     1262
##    y0[10]  28.06  28.07 4.84 4.41  20.36  35.22 1.00     1451     1271
##    y0[11]  26.64  26.67 4.59 4.15  19.08  34.33 1.00     1850     1891
##    y0[12]  26.98  27.14 4.88 4.24  19.14  34.82 1.00     1842     1593
##    y0[13]  27.82  27.86 4.67 4.16  20.50  35.36 1.00     1415     1186
##    m1[1]   17.46  17.43 3.05 2.64  12.98  21.70 1.03      264      316
##    m1[2]   18.83  18.81 2.55 2.16  14.97  22.38 1.03      270      326
##    m1[3]   20.20  20.21 2.07 1.80  17.00  23.18 1.02      290      307
##    m1[4]   21.58  21.60 1.65 1.44  18.97  24.06 1.02      341      329
##    m1[5]   22.95  22.95 1.34 1.20  20.77  25.02 1.01      541      614
##    m1[6]   24.32  24.34 1.22 1.10  22.35  26.21 1.00     1216      978
##    m1[7]   25.69  25.74 1.35 1.24  23.50  27.86 1.00     1674      928
##    m1[8]   27.07  27.12 1.67 1.45  24.32  29.67 1.00      809      710
##    m1[9]   28.44  28.48 2.09 1.83  25.01  31.68 1.01      554      575
##    m1[10]  29.81  29.83 2.57 2.28  25.71  33.77 1.01      457      582
##    y1[1]   17.46  17.46 5.37 4.65   9.06  26.00 1.01      629      771
##    y1[2]   18.74  18.77 4.99 4.50  10.96  26.75 1.01      823      804
##    y1[3]   20.08  20.18 4.81 4.51  12.50  27.73 1.00     1135     1265
##    y1[4]   21.58  21.63 4.77 4.14  13.90  29.31 1.01     1047     1630
##    y1[5]   23.08  23.05 4.55 4.37  16.10  30.64 1.00     1449     1428
##    y1[6]   24.30  24.21 4.31 4.06  17.62  31.16 1.00     1820     1738
##    y1[7]   25.69  25.58 4.50 3.98  18.48  33.11 1.00     2025     1510
##    y1[8]   27.08  27.17 4.53 4.16  19.58  34.17 1.00     1703     1606
##    y1[9]   28.30  28.27 4.70 4.27  20.56  35.96 1.00     1780     1776
##    y1[10]  29.92  29.94 5.15 4.59  21.62  38.20 1.01     1044     1167

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -99.1728 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -28.7063   0.000276867   2.73144e-05           1           1       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -28.71
##    b0        11.13
##    b1         2.46
##    s          4.08
##    m0[1]     21.57
##    m0[2]     20.76
##    m0[3]     14.02
##    m0[4]     22.23
##    m0[5]     22.35
##    m0[6]     15.07
##    m0[7]     26.61
##    m0[8]     18.96
##    m0[9]     29.69
##    m0[10]    29.69
##    m0[11]    27.69
##    m0[12]    27.85
##    m0[13]    29.32
##    y0[1]     26.99
##    y0[2]     22.04
##    y0[3]     16.42
##    y0[4]     21.22
##    y0[5]     20.17
##    y0[6]     17.13
##    y0[7]     27.84
##    y0[8]     20.35
##    y0[9]     23.36
##    y0[10]    32.12
##    y0[11]    26.68
##    y0[12]    25.98
##    y0[13]    31.74
##    m1[1]     12.31
##    m1[2]     14.56
##    m1[3]     16.81
##    m1[4]     19.06
##    m1[5]     21.31
##    m1[6]     23.56
##    m1[7]     25.81
##    m1[8]     28.06
##    m1[9]     30.31
##    m1[10]    32.56
##    y1[1]     11.28
##    y1[2]     10.30
##    y1[3]     17.24
##    y1[4]     20.66
##    y1[5]     18.76
##    y1[6]     14.29
##    y1[7]     25.87
##    y1[8]     26.28
##    y1[9]     23.49
##    y1[10]    32.09
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -29.03 -28.64 1.44 1.14 -32.08 -27.49 1.01      424      553
##    b0      10.10  10.41 3.10 2.76   4.50  14.27 1.02      237      253
##    b1       2.65   2.58 0.57 0.51   1.87   3.73 1.01      293      258
##    s        5.16   4.87 1.41 1.22   3.43   7.76 1.00      813      932
##    m0[1]   21.36  21.45 1.37 1.22  18.97  23.43 1.00      504      588
##    m0[2]   20.48  20.59 1.44 1.28  17.89  22.59 1.00      405      457
##    m0[3]   13.22  13.50 2.51 2.22   8.70  16.62 1.02      239      260
##    m0[4]   22.07  22.15 1.33 1.22  19.77  24.09 1.00      629      692
##    m0[5]   22.20  22.26 1.33 1.21  19.90  24.20 1.00      655      719
##    m0[6]   14.35  14.59 2.31 2.07  10.16  17.53 1.01      242      274
##    m0[7]   26.80  26.78 1.51 1.38  24.39  29.37 1.01     1473     1131
##    m0[8]   18.55  18.72 1.65 1.45  15.49  20.92 1.01      300      331
##    m0[9]   30.12  30.03 1.96 1.83  27.11  33.54 1.01      905      909
##    m0[10]  30.12  30.02 1.96 1.83  27.11  33.54 1.01      906      909
##    m0[11]  27.96  27.92 1.65 1.53  25.34  30.79 1.01     1275      990
##    m0[12]  28.14  28.09 1.67 1.55  25.46  31.03 1.01     1240      966
##    m0[13]  29.72  29.63 1.90 1.79  26.80  33.05 1.01      968      875
##    y0[1]   21.42  21.59 5.43 5.00  12.48  30.22 1.00     2006     1586
##    y0[2]   20.28  20.46 5.62 5.29  10.93  29.15 1.00     1928     1778
##    y0[3]   13.20  13.34 5.96 5.51   3.32  22.70 1.00     1250     1329
##    y0[4]   21.94  22.05 5.35 4.91  13.43  30.42 1.00     1822     1973
##    y0[5]   22.48  22.56 5.77 5.08  13.32  31.31 1.00     2020     1431
##    y0[6]   14.52  14.80 5.89 5.35   4.23  23.33 1.00     1109     1231
##    y0[7]   26.77  26.72 5.46 5.02  17.84  35.73 1.00     1983     1913
##    y0[8]   18.43  18.53 5.60 5.01   9.24  27.24 1.00     1394     1568
##    y0[9]   30.10  30.03 5.62 5.01  21.27  39.47 1.00     1813     1895
##    y0[10]  30.07  29.96 5.84 5.19  20.34  39.64 1.00     1797     1754
##    y0[11]  27.81  27.72 5.77 5.09  18.77  37.37 1.00     1792     1750
##    y0[12]  28.24  28.27 5.57 5.18  19.23  37.26 1.00     2034     1698
##    y0[13]  29.70  29.54 5.67 5.26  20.70  39.36 1.00     1886     1242
##    m1[1]   11.37  11.67 2.86 2.53   6.23  15.22 1.02      237      261
##    m1[2]   13.80  14.06 2.41 2.14   9.46  17.08 1.01      240      273
##    m1[3]   16.23  16.45 1.99 1.79  12.52  19.01 1.01      255      269
##    m1[4]   18.65  18.83 1.64 1.43  15.62  21.02 1.01      303      342
##    m1[5]   21.08  21.17 1.39 1.23  18.64  23.14 1.00      465      588
##    m1[6]   23.51  23.51 1.31 1.21  21.35  25.46 1.01     1032      990
##    m1[7]   25.94  25.91 1.43 1.29  23.67  28.34 1.01     1509     1139
##    m1[8]   28.36  28.31 1.70 1.59  25.67  31.33 1.01     1197      976
##    m1[9]   30.79  30.67 2.07 1.91  27.61  34.46 1.01      819      904
##    m1[10]  33.22  33.02 2.50 2.34  29.51  37.76 1.01      629      601
##    y1[1]   11.38  11.77 6.02 5.29   0.89  20.60 1.00      864     1156
##    y1[2]   13.97  14.23 6.02 5.51   3.79  23.28 1.00     1089     1142
##    y1[3]   16.17  16.33 5.81 5.12   6.35  25.35 1.00     1174     1172
##    y1[4]   18.58  18.70 5.73 5.18   9.36  27.49 1.00     1207     1298
##    y1[5]   21.20  21.20 5.57 4.99  11.90  30.27 1.00     1685     1654
##    y1[6]   23.52  23.55 5.51 5.01  14.48  32.31 1.00     1902     1961
##    y1[7]   25.94  25.82 5.61 5.26  16.85  35.08 1.00     1895     1638
##    y1[8]   28.34  28.30 5.56 5.08  19.67  37.42 1.00     1760     1434
##    y1[9]   30.83  30.67 5.62 5.16  21.97  39.99 1.00     1297     1379
##    y1[10]  33.20  33.23 5.86 5.24  23.78  42.71 1.00     1725     1782

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -16.7538 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -13.7528    0.00030773    2.0269e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -13.75
##      p        0.78
##      se       0.54
##      sp       0.43
##      ppv      0.77
##      npv      0.21
system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,ch=F)
})
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -19.54 -19.20 1.29 1.01 -22.07 -18.13 1.00      718      912
##      p      0.51   0.50 0.29 0.37   0.05   0.95 1.01      609      458
##      se     0.53   0.53 0.12 0.12   0.34   0.72 1.00     3217     1505
##      sp     0.45   0.44 0.16 0.17   0.20   0.71 1.00     2121     1464
##      ppv    0.50   0.50 0.30 0.39   0.05   0.95 1.01      632      479
##      npv    0.48   0.47 0.29 0.38   0.05   0.95 1.01      633      539

##   ユーザ システム     経過 
##    0.558    0.121    0.699

appendix: bimodal distribution, mixed normal distribution

stan

data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real<lower=0, upper=1> theta; //ratio of mix
  real mu1;
  real mu2;
  real<lower=0> sigma1;
  real<lower=0> sigma2;
}
model {
  for (n in 1:N) {
    target += log_mix(theta,
                      normal_lpdf(y[n] | mu1, sigma1),
                      normal_lpdf(y[n] | mu2, sigma2));
  }
}

EM algorithm

y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -239 100  6 -506 -507
## 
## Clustering table:
##  1  2  3 
## 28 44 28
rst$parameters
## $pro
## [1] 0.281 0.438 0.281
## 
## $mean
##      1      2      3 
## -5.214 -0.047  5.029 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.842
## 
## 
## $Vinv
## NULL
plot(rst)